home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / TMTJ.CPP < prev    next >
C/C++ Source or Header  |  1995-01-11  |  5KB  |  155 lines
  1.  
  2. //#define WANT_STREAM
  3.  
  4. #include "include.h"
  5.  
  6. #include "newmatap.h"
  7. //#include "newmatio.h"
  8.  
  9. void Print(const Matrix& X);
  10. void Print(const UpperTriangularMatrix& X);
  11. void Print(const DiagonalMatrix& X);
  12. void Print(const SymmetricMatrix& X);
  13. void Print(const LowerTriangularMatrix& X);
  14.  
  15. void Clean(Matrix&, Real);
  16.  
  17.  
  18.  
  19.  
  20. void trymatj()
  21. {
  22.    Tracer et("Nineteenth test of Matrix package");
  23.    Exception::PrintTrace(TRUE);
  24.    // testing elementwise products
  25.  
  26.    {
  27.       Tracer et1("Stage 1");
  28.       Matrix A(13,7), B(13,7), C(13,7);
  29.       int i,j;
  30.       for (i=1;i<=13;i++) for (j=1; j<=7; j++)
  31.       {
  32.           Real a = (i+j*j)/2, b = (i*j-i/4);
  33.           A(i,j)=a; B(i,j)=b; C(i,j)=a*b;
  34.       }
  35.       // Where complete matrix routine can be used
  36.       Matrix X = SP(A,B)-C; Print(X);
  37.       X = SP(A,B+1.0)-A-C; Print(X);
  38.       X = SP(A-1,B)+B-C; Print(X);
  39.       X = SP(A-1,B+1)+B-A-C+1; Print(X);
  40.       // Where row-wise routine will be used
  41.       A = A.Rows(7,13); B = B.Rows(7,13); C = C.Rows(7,13);
  42.       LowerTriangularMatrix LTA; LTA << A;
  43.       UpperTriangularMatrix UTB; UTB << B;
  44.       DiagonalMatrix DC; DC << C;
  45.       X = SP(LTA,UTB) - DC; Print(X);
  46.       X = SP(LTA*2,UTB) - DC*2; Print(X);
  47.       X = SP(LTA, UTB /2) - DC/2; Print(X);
  48.       X = SP(LTA/2, UTB*2) - DC; Print(X);
  49.       DiagonalMatrix DX;
  50.       DX << SP(A,B); DX << (DX-C); Print(DX);
  51.       DX << SP(A*4,B); DX << (DX-C*4); Print(DX);
  52.       DX << SP(A,B*2); DX << (DX-C*2); Print(DX);
  53.       DX << SP(A/4,B/4); DX << (DX-C/16); Print(DX);
  54.       LowerTriangularMatrix LX;
  55.       LX = SP(LTA,B); LX << (LX-C); Print(LX);
  56.       LX = SP(LTA*3,B); LX << (LX-C*3); Print(LX);
  57.       LX = SP(LTA,B*5); LX << (LX-C*5); Print(LX);
  58.       LX = SP(-LTA,-B); LX << (LX-C); Print(LX);
  59.    }
  60.    {
  61.       // Symmetric Matrices
  62.       Tracer et1("Stage 2");
  63.       SymmetricMatrix A(25), B(25), C(25);
  64.       int i,j;
  65.       for (i=1;i<=25;i++) for (j=i;j<=25;j++)
  66.       {
  67.          Real a = i*j +i - j + 3;
  68.          Real b = i * i + j;
  69.          A(i,j)=a; B(i,j)=b; C(i,j)=a*b;
  70.       }
  71.       UpperTriangularMatrix UT;
  72.       UT << SP(A,B); UT << (UT - C); Print(UT);
  73.       Matrix MA = A, X;
  74.       X = SP(MA,B)-C; Print(X);
  75.       X = SP(A,B)-C; Print(X);
  76.       SymmetricBandMatrix BA(25,5), BB(25,5), BC(25,5);
  77.       BA.Inject(A); BB.Inject(B); BC.Inject(C);
  78.       X = SP(BA,BB)-BC; Print(X);
  79.       X = SP(BA*7,BB)-BC*7; Print(X);
  80.       X = SP(BA,BB/8)-BC/8; Print(X);
  81.       X = SP(BA*16,BB/16)-BC; Print(X);
  82.       X = SP(BA,BB); X=X-BC; Print(X);
  83.       X = SP(BA*2, BB/2)-BC; Print(X);
  84.       X = SP(BA, BB/2)-BC/2; Print(X);
  85.       X = SP(BA*2, BB)-BC*2; Print(X);
  86.    }
  87.    {
  88.       // Band matrices
  89.       Tracer et1("Stage 3");
  90.       Matrix A(19,19), B(19,19), C(19,19);
  91.       int i,j;
  92.       for (i=1;i<=19;i++) for (j=1;j<=19;j++)
  93.       {
  94.          Real a = i*j +i - 1.5*j + 3;
  95.          Real b = i * i + j;
  96.          A(i,j)=a; B(i,j)=b; C(i,j)=a*b;
  97.       }
  98.       BandMatrix BA(19,10,7), BB(19,8,15), BC(19,8,7);
  99.       BA.Inject(A); BB.Inject(B); BC.Inject(C);
  100.       Matrix X; BandMatrix BX; ColumnVector BW(2);
  101.       X = SP(BA,BB); X=X-BC; Print(X);
  102.       X = SP(BA/8,BB); X=X-BC/8; Print(X);
  103.       X = SP(BA,BB*17); X=X-BC*17; Print(X);
  104.       X = SP(BA/4,BB*7); X=X-BC*7/4; Print(X);
  105.       X = SP(BA,BB)-BC; Print(X);
  106.       X = SP(BA/8,BB)-BC/8; Print(X);
  107.       X = SP(BA,BB*17)-BC*17; Print(X);
  108.       X = SP(BA/4,BB*7)-BC*7/4; Print(X);
  109.       BX = SP(BA,BB);
  110.       BW(1)=BX.upper-7; BW(2)=BX.lower-8; Print(BW);
  111.  
  112.       BA.ReDimension(19,7,10); BB.ReDimension(19,15,8);
  113.       BC.ReDimension(19,7,8);
  114.       BA.Inject(A); BB.Inject(B); BC.Inject(C);
  115.  
  116.       X = SP(BA,BB); X=X-BC; Print(X);
  117.       X = SP(BA/8,BB); X=X-BC/8; Print(X);
  118.       X = SP(BA,BB*17); X=X-BC*17; Print(X);
  119.       X = SP(BA/4,BB*7); X=X-BC*7/4; Print(X);
  120.       X = SP(BA,BB)-BC; Print(X);
  121.       X = SP(BA/8,BB)-BC/8; Print(X);
  122.       X = SP(BA,BB*17)-BC*17; Print(X);
  123.       X = SP(BA/4,BB*7)-BC*7/4; Print(X);
  124.       BX = SP(BA,BB);
  125.       BW(1)=BX.upper-8; BW(2)=BX.lower-7; Print(BW);
  126.    }
  127.    {
  128.       // SymmetricBandMatrices
  129.       Tracer et1("Stage 4");
  130.       Matrix A(7,7), B(7,7);
  131.       int i,j;
  132.       for (i=1;i<=7;i++) for (j=1;j<=7;j++)
  133.       {
  134.          Real a = i*j +i - 1.5*j + 3;
  135.          Real b = i * i + j;
  136.          A(i,j)=a; B(i,j)=b;
  137.       }
  138.       BandMatrix BA(7,2,4), BB(7,3,1), BC(7,2,1);
  139.       BA.Inject(A);
  140.       SymmetricBandMatrix SB(7,3);
  141.       SymmetricMatrix S; S << (B+B.t());
  142.       SB.Inject(S); A = BA; S = SB;
  143.       Matrix X;  
  144.       X = SP(BA,SB); X=X-SP(A,S); Print(X);
  145.       X = SP(BA*2,SB); X=X-SP(A,S*2); Print(X);
  146.       X = SP(BA,SB/4); X=X-SP(A/4,S); Print(X);
  147.       X = SP(BA*4,SB/4); X=X-SP(A,S); Print(X);
  148.       X = SP(BA,SB)-SP(A,S); Print(X);
  149.       X = SP(BA*2,SB)-SP(A,S*2); Print(X);
  150.       X = SP(BA,SB/4)-SP(A/4,S); Print(X);
  151.       X = SP(BA*4,SB/4)-SP(A,S); Print(X);
  152.    }
  153.  
  154. }
  155.